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ABSTRACT 

Three-dimensional turbulent jets with rectangular cross-section are simulated with a finite-dif- 
ference numerical method. The full Navier-Stokes equations are solved at low Reynolds numbers, 
whereas at the high Reynolds numbers filtered forms of the equations are solved along with a sub- 
grid scale model to approximate effects of the unresolved scales. A 2-N storage, third-order 
Runge-Kutta scheme is used for temporal discretization and a fourth-order compact scheme is 
used for spatial discretization. Computations are performed for different inlet conditions which 
represent different types of jet forcing. The phenomenon of axis- switching is observed, and it is 
confirmed that this is based on self-induction of the vorticity field. Budgets of the mean stream- 
wise velocity show that convection is balanced by gradients of the Reynolds stresses and the pres- 
sure. 
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INTRODUCTION 


Turbulent jets are present in many physical processes and technological applications. Turbu- 
lent jets can be found in combustors where the fuel and oxidizer are introduced as co-flowing jets. 
The efficiency of such a process is largely determined by the mixing of the jets. Recently, jet air- 
craft noise has received much attention due to plans for a high-speed civil transport vehicle. A 
critical issue for the project’s success is reducing jet noise to acceptable levels near populated 
areas. The belief is that acoustic patterns can be altered by manipulating the large scale structures 
in turbulent jet flows through external forcing. Non-circular jets can also be used to enhance the 
mixing of hot jet gases with the surroundings and thus avoid aircraft detection. In industrial appli- 
cations, efficient mixing is required to dilute pollutants issuing from smokestacks with the ambi- 
ent air to minimize its harmful effects or to promote fuel-oxidant mixing in combustion chambers. 

Experiments (Tsuchiya et al. 1986, 1989, Quinn 1989, 1994) have shown that three-dimen- 
sional (3-D) jets can be used to enhance mixing and entrainment rates in comparison to nominally 
two-dimensional (2-D) jets. A fundamental understanding of the dynamics of complex, turbulent 
jets is required for their prediction and control. The present study is concerned with the under- 
standing of the role of vorticity in the spatial evolution of incompressible 3-D jets in the near to 
medium field and the effects of external forcing thereupon. Previous numerical simulations of 
such flows have been performed by Grinstein (1995), with emphasis on the vorticity dynamics in 
the near field (less than 5 diameters), and Miller et al. (1995), who studied co-flowing jets at a low 
Reynolds number of 800. In the latter, examination of intantaneous vorticity fields show symme- 
tries which are not indicative of turbulent flow. This work extends those studies into the medium 
field, beyond the extent of the potential core, and into truly turbulent flows. 

MATHEMATICAL FORMULATION 

The partial differential equations governing the incompressible jet fluid flow are the Navier- 
Stokes equations which can be written in Cartesian tensor form, for dimensionless variables as: 
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where, Uj are the Cartesian velocity components in the Cartesian coordinate directions Xj, p is the 
pressure and Re D * is the Reynolds number based on the equivalent diameter. These equations 
must be solved in conjuction with the continuity equation: 


dll: 


dx 


- = 0 


which expresses the divergence-free velocity condition. 


( 2 ) 


At higher Reynolds numbers, all scales present in the flow cannot be resolved on computa- 
tional grids that present resources would allow, so a large eddy simulation (LES) must be utilized. 
The application of a grid filter, G, to equations (1) and (2) results in the filtered equations of 
motion, in which the overbar indicates a filtered variable: 


and 
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where u- represents components of the filtered velocity which is resolved in the computation, and 
Xjj = p UjUj - pUjlij is the subgrid scale (SGS) Reynolds stress which must be modeled in terms 
of the resolved velocity field. The Smagorinsky eddy-viscosity model is utilized in the present 
study to approximate the SGS Reynolds stress as: 


^ 3 




where A is a length scale associated with the grid size, Su = 


f du i diAj] 
y dXj dx , J 


(5) 

is the resolved strain 


rate tensor, and |sJ = j2SjjSjj . The model coefficient, C 5 , is set to the constant value of 0.01 . 
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The Navier-Stokes equations (1) and the filtered equations (3) are discretized temporally with 
explicit Runge-Kutta (RK) schemes and spatially with implicit compact finite difference schemes. 
The discretized equation has the form: 


n + 1 n j M i r y t M o i 
M. = Uj +b At[H i -0 P ] 


with 
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i.e., the sum of convection and diffusion (CD) terms. 


( 6 ) 


In Eq.(6), n represents the time step and M stands for the M th stage of the RK scheme, with the 
corresponding coefficient b M . 5 X and 5 rx . are compact first and second derivative operators, 
respectively. In the present study low-storage RK schemes are utilized. The low-storage require- 
ment is accomplished by continuously overwriting the storage location for the time derivatives 
and unknown variables at each sub-stage: 



M S, M 
< — Cl H : 


M + 1 M ,M . -M 

u < — u, + b AtH : 


(7) 

(8) 


where H™ = - [dp M /dxj\ and the notation is used to indicate that the storage 

locations, H f 1 , ii* are overwritten by, (if , + 1 , respectively. The coefficients a M and b M 

for the low-storage, third-order scheme (Williamson, 1980; Lowery and Reynolds 1986) are. 


Table 1 : Coefficients of the third-order RK scheme 


M 

a M 

h M 

1 

0 

0.500 

2 

-0.68301270 

0.91068360 

3 

-1.33333333 

0.36602540 
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Compact finite-difference schemes (Leie 1992) used to approximate derivatives which appear 
in Eqs. (1-5) involve implicit treatment of discrete derivatives and explicit treatment of discrete 
functions. Only tri-diagonal schemes are utilized in this study. For example: 


«<!>', - i + <t>', + + , = 2I -(<►/ + 1 - - 1 > ( 9 ) 

where 0’- represents the first derivative of the generic variable 0 with respect to x-, and a, are 
the coefficients of the compact scheme which determine the accuracy. For the fourth-order 
scheme: a = 1/4 ; a = 3/2 . The LHS of Eq. (9) contains the unknown derivatives at grid 
points i and i ± 1 while the RHS contains the known functional values 0 • at the grid points i ± 1 . 
Similarly, the second derivative terms present in the viscous terms of the momentum equation and 
the Laplacian operator of the Poisson equation for pressure are approximated using fourth-order 
compact finite differences: 

a0" ( _| +0” / + a0” l+1 = — ^(0 I+1 +20 ( + 0._,) (10) 

(Ax ,.) 2 

where 0", represents the second derivative of the generic variable 0, with respect to x, , and a , a, 
b are the coefficients of the compact scheme. For the fourth-order scheme: a = 1 / 10 ; 
a = 6/5 . The tri-diagonal system of algebraic equations are solved using the Thomas algorithm. 
At the boundary third-order, one-sided differences are used to close the system of equations aris- 
ing from the first or second derivative schemes. 

A major requirement in the computation of incompressible flow is the satisfaction of the 
divergence-free velocity condition. This is made particularly difficult by the absence of an evolu- 
tion equation for pressure. This condition must be satisfied indirectly through the solution of a 
Poisson equation derived by taking the divergence of the momentum equations: 


V 2 P M = 8. 
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The discretization of the Laplacian operator V results in: 
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( 12 ) 




M 
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It is necessary to solve this equation completely at every substage of the RK scheme. Although 
the computational matrix is sparse, direct methods of solutions are not feasible. An iterative multi- 
grid method was found to be most efficient, with reasonably good convergence rates. Neverthe- 
less, the solution of the Poisson equation accounts for about 50% of the total computational effort. 


In order to solve a well posed problem, the boundary and initial conditions for the jet simula- 
tions are defined. The ellipticity in the spatial terms of the governing equations requires that 
boundary conditions be defined on all boundaries. In the laboratory, jet flows are commonly gen- 
erated by the use of a fan which forces fluid along an enclosed duct of nozzle. The jet leaves the 
exit plane of the nozzle where it interacts with the ambient fluid. Prior to exit, the jet can be con- 
sidered as a relatively uniform freestream and a curved boundary layer at the walls of the nozzle. 
A short distance downstream of the nozzle exit, the boundary layer is smoothened so that the 
streamwise velocity can be modeled using the hyperbolic tangent (tanh) function. The tanh func- 
tion enables the streamwise velocity to transition in the radial direction from the uniform velocity 
at the core of the jet to that of the ambient fluid in the freestream. The inflow boundary of the 
computational domain is placed at a short distance downstream of the nozzle exit which is not 
actually included in the jet simulations. 

The mean or time-averaged streamwise velocity component at the inflow boundary is given 
by: 


- ( U H -U .) ( r \ 

«,( 0,X 2 ,X 3 ) = U c + tanh^— j (13) 

where U c = (U H + U L )/ 2 is the convective velocity and U H , U L , represent the velocities of the 
jet core, and ambient fluid, respectively. The variable r, represents the minimum directed distance 
from the point, ( 0 , x 2 , x 3 ) to the line of constant convective velocity of the boundary layer. The 
momentum thickness of the boundary layer at the inflow plane, 0^ is used to normalize the 
directed distance, r. If the point (0, x 2 , x 3 ) is “outside” the boundary layer, r is defined to be nega- 
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tive, otherwise r is defined to be positive in the point lies on the inside of the boundary layer. 
Equation (13) produces a constant thickness boundary layer if the momentum thickness, 0 O , is 
constant at all azimuthal positions along the boundary layer. Non constant thickness boundary 
layers can be generated by specifying the desired variation of B 0 along the perimeter of the 
boundary layer. The mean major and minor direction velocity components at the inflow plane, 
u 2 (0, x 2 , x 3 ) and u 3 (0, x 2 , x 3 ) are specified from experiment. 

To simulate instabilities in the boundary layer, a time dependent forcing function of low inten- 
sity is added to the mean velocity components at the inflow boundary to promote unsteady 
motion. At higher Reynolds numbers or large computational lengths small round-off errors would 
grow to produce unsteady motion of the unstable shear layers thus obviating the need for forcing 
functions. Two classes of perturbations are used in the current study; (/) sinusoidal perturbations 
with frequency and amplitudes related to the most unstable modes found from viscous stability 
analysis (Wilson and Demuren 1996); (ii) perturbations having an experimentally measured 
velocity spectrum and transverse root mean square ( rms ) value. Perturbations which have a broad 
spectrum resembling that of fully-developed, mostly random, 3-D turbulence were generated 
(Wilson 1993). The perturbations are typical of those found in the experimental jets originating 
from contoured nozzles. The power spectrum and rms values are taken from experiment. Because 
phase information is not included in the power spectrum, a random phase relationship for the 
modes comprising the spectrum was assumed. The random inlet boundary conditions are the spa- 
tial analog to the random perturbations generated for initial conditions in temporal simulations. 
Class (i) represents jets issuing from a nozzle with laminar boundary layers, whereas class (ii) 
represents jets issuing from a nozzle with turbulent boundary layers. 

MODEL PROBLEMS 

Spatial simulations of rectangular jet flows are performed in this study in which a fixed region 
of the flow is computed and disturbances grow in the streamwise direction. This can be contrasted 
with a temporal simulation in which a small region of the flow is followed in time and the domain 
moves in the streamwise direction. Spatial simulations are closer to reality, but are computation- 
ally more demanding. As a result of the spatial reference frame, initial conditions are of minor 
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importance because they are quickly convected out of the domain and the dynamics of the jet flow 
are determined by the forcing functions applied at the inflow plane, until the inherent instabilities 
of the flow take over or the onset of transition to turbulent flow. Simulations are started on coarser 
grids with the velocities specified at the inflow plane used to initialize the velocity field in the inte- 
rior at t = 0. After several flow through times (the time required for a fluid particle to convect from 
the inflow to the outflow plane at the jet core velocity, U e ) the initial conditions are “washed” 
from the domain. Simulations on finer grids are commenced with initial fields obtained by prolon- 
gating the coarser grid results using a standard, 2nd-order accurate interpolation formula. 

The rectangular jet of the present study has a nominal aspect ratio of 2:1. This corresponds to 
the configuration of some jet issuing from contoured nozzle studied by Quinn (private communi- 
cation). Figure 1(c) illustrates the shape of the jet cross-section at exit. Cases with both low and 
high Reynolds numbers are considered. In the former, the flow is fully resolved and no sub-grid 
scale model is required, but in the latter, filtered equations are solved and an SGS model is uti- 
lized. Reynolds numbers based on the core velocity and the equivalent diameter are 750 and 
75000, respectively. The four cases presented in this paper correspond to the two Reynolds num- 
bers and the different types of inlet forcing functions, as explained previously, namely; (i) sinuso- 
idal perturbations with the fundamental mode at an rms of 3% of mean jet velocity, and Re = 750; 
(ii) sinusoidal perturbation with the fundamental and first subharmonic modes, each at an rms of 
1.5%, and Re = 750; (iii) same perturbations as in (ii), but at Re = 75000; (iv) random broad-mode 
perturbations with an rms of 3% of mean jet velocity, and Re = 75000. Case (i) was simulated in a 
domain with size (10 x 5 x 5) on a uniform computational grid of (80 x 65 x 65), whereas all other 
cases were simulated on a (10 x 10 x 10) domain with computational grids of (100 x 129 x 129). 
The computational grids correspond to 5 points per wavelength (PPW) of the fundamental shear 
layer instability, in the transverse direction, and 10 PPW in the streamwise direction. Analysis of 
the fourth-order compact scheme in use here, in several benchmark problems (Wilson 1996), has 
shown that roughly 6 PPW are sufficient for accurate simulation. The compact scheme is inher- 
ently non-diffusive and non-dissipative, and with adequate resolution dispersive errors are mini- 
mized. 
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RESULTS AND DISCUSSION 

Figures 1 and 2 show computed instantaneous vorticity fields for cases (i) and (ii), respec- 
tively. In case (i), the flow remains symmetrical about the major and minor axes, just like the inlet 
mean flow and perturbation velocities, hence the flow has remained laminar and transition to tur- 
bulence did not occur in the 10 diameters computed. The phenomenon of axis switching is 
observed to be taking place near the end of the domain. The source of this is seen to be clearly in 
the complex vorticity field. It is most likely due to self-induction, since there is no streamwise 
vorticity at the inlet plane. It appears that the effect of the induced streamwise vorticity field is to 
pump fluid from the major axis to the minor axis. Now, since the flow remains laminar, this pro- 
cess is inviscid and is unrelated to turbulence-generated secondary flow as found in streamwise 
corners, the so-called secondary flow of Prandtl’s second kind, (e.g., Demuren and Rodi 1984). 
The role of viscosity is to promote the eventual decay of this secondary motion. In case (ii), the 
vorticity field is even more complex, indicating the effects of the sub-harmonic component in 
inducing multiple vortex pairings. Near the end of the computed domain, some asymmetry can be 
discerned, indicating the beginning of the transition process. In this case, the phenomenon of axis- 
switching was not observed. It is therefore quite clear that the modes of instability waves present 
in the emerging boundary layer have strong effects on the evolution of the jet. 

Time-averaged jet half-widths are shown for the two high Reynolds number cases, (iii) and (v) 
in Fig. 3. We see that in case (iii), with discrete two-mode forcing, no axis-switching occurred up 
to ten diameters downstream. This result is in agreement with that of case (ii), and is further con- 
firmation of the thesis that the phenomenon of axis-switching is not controlled by the turbulence 
structure. On the other hand, in case (iv) axis-switching occurred at about 7 diameters. This inlet 
condition, with a broad-mode spectrum of perturbations is indicative of a jet evolving from a noz- 
zle with turbulent boundary layers. Experimental data from similar conditions (Tsuchiya and 
Horikoshi 1986) suggest that axis-switching was observed at between 6 to 8 diameters. However, 
the presence of axis-switching should not be equated directly with increased entrainment of ambi- 
ent fluid which is desired in several engineering applications. Integration of the entrained flow 
shows that more ambient fluid was entrained in case (iii) than in case (iv). 
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All subsequent results focus on case (iv), which is closest to experiment of unforced jets. Fig- 
ures 4 and 5 show instantaneous fields of the streamwise velocity, pressure, total vorticity, and the 
streamwise vorticity in the major and minor axis planes, respectively, after nine flow through 
times. Displayed numbers correspond to peak values for each variable, and negative values are 
drawn with dashed lines. The potential core of the jet extends to about 4 diameters, and significant 
streamwise vorticity commences from this location downstream. The latter quickly reached peak 
values by about 6 to 7 diameters, and subsequently decayed in strength further downstream. The 
location of axis cross-over corresponds roughly to the location of peak streamwise vorticity. Fur- 
thermore, although there are some pockets of positive instantaneous pressure, it is largely nega- 
tive, thus fueling the entrainment of ambient fluid. 

Figure 6 shows the balance of terms in the streamwise momentum equation, time-averaged 
over 1000 time steps and approximately 10 flow through times. We see that, by an overwhelming 
margin, the spatial convection is balanced by gradients of the resolved Reynolds stresses. By com- 
parison, gradients of the SGS stresses are negligible, which would indicate the reliability of the 
LES approach for this problem even at the chosen Reynolds number of 75000. The third most 
important term is the pressure gradient. In fact, near the edges, it is the dominant term which bal- 
ances spatial convection, i.e., it acts as the source for the entrainment of ambient fluid. The sym- 
bols, in the figure, represent the lack of balance or drift from steady state, which is indicative of 
the adequacy of the sample size used in the time-averaging process. With no homogeneous direc- 
tion in the present flow, the sample size seems to be adequate. Further time-averaged quantities of 
first and second moments, as well as budgets of the respective equations can be found in Wilson 
(1996). 

CONCLUDING REMARKS 

Numerical simulation of jets with rectangular cross-section have been performed using a 
third-order Runge-Kutta scheme for temporal integration and a fourth-order compact scheme for 
spatial discretization. Cases at low Reynolds numbers are direct simulations and those at high 
Reynolds number are large eddy simulations. 
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The phenomenon of axis-switching is observed to be dependent on instability waves present 
in the inlet boundary layers, and could be induced in both laminar and turbulent jets. It is con- 
firmed that this is based on self-induction of the vorticity field. The presence of a discrete sub-har- 
monic mode led to its suppression, but entrainment of ambient fluid was still enhanced. 

Budget of terms in the mean streamwise momentum equation showed that spatial convection 
is balanced mostly by gradients of the Reynolds stresses, with pressure gradients making a minor 
contribution. However, near the edges, ambient fluid entrainment is fueled by the pressure gradi- 
ents. Terms involving the SGS stresses and viscous diffusion are negligible. 
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Figure 1. Contours of vorticity magnitude, (a) - (e) and streamwise vorticity, (f) - (h) for case (i) at 
t = 2 flow through times, for fundamental forcing function, (a) minor axis plane, z/D e = 
0, (b) major axis plane, y/D e = 0, (c) cross flow plane, x/D e = 0, (d) x/D e = 5, (e) x/D e = 
10, (f) x/D e = 2.5, (g) x/D e = 5, and (h) xfD e = 10. 
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Figure 2. Contours of vorticity magnitude, (a) - (e) and streamwise vorticity, (f) - (h) for case (ii) 
at t = 2 flow through times, for fundamental and first subharmonic forcing function. (a) 
minor axis plane, z/D e - 0, (b) major axis plane, y/D e = 0, (c) cross flow plane, x/D e = 0, 
(d) x/D e = 5, (e) x/D e = 10, (f) x/D e = 2.5, (g) x/D e = 5, and (h) x/D e = 10. 








Figure 3. Jet half-widths for (a) case (iii), (b) case (iv) (solid - major axis plane, dashed - minor 
axis plane) . 
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Figure 4. Instantaneous contours of dimensionless (a) streamwise velocity, (b) pressure, (c) total 
vorticity, (d) streamwise vorticity along major axis plane at t = 9 flow through times, for 
case (iv). (Numbers represent peak values). 
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Figure 5. Instantaneous contours of dimensionless (a) streamwise velocity, (b) pressure, (c) total 
vorticity, (d) streamwise vorticity along minor axis plane, at t = 9 flow through times, for 
case (iv). (Numbers represent peak values). 
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